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Based on a simple geometrical approach, we analyze the evolution of the Kozai-Lidov mechanism 
for stars around shrinking massive black hole binaries on circular orbits. We find that, due to 
a peculiar bifurcation pattern induced by the Newtonian potential of stellar clusters, the orbit 
of stars could become highly eccentric. This transition occurs abruptly for stars with small initial 
eccentricities. Our approach would be also useful for studying the Kozai-Lidov mechanism in various 
astrophysical contexts. 


I. INTRODUCTION 

Coalescence of a massive black hole (MBH) binary is 
a spectacular event and a huge amount of energy is con¬ 
verted into gravitational waves (GWs). Depending on 
the masses of the binary, the planned space interferom¬ 
eters such as eLISA [I| can detect the strong GW sig¬ 
nal, even from a very high redshift. Meanwhile, associ¬ 
ated with shrinkage of MBH binaries in galactic nuclei, 
the infall rate of stars (more precisely, stellar-mass ob¬ 
jects) to the MBHs could be enhanced significantly 0. 
These stars would emit GWs or electromagnetic waves 
that encode precious information for astrophysics, cos¬ 
mology and fundamental physics {e.g. strong gravity) 

i- 

For the infall into the MBHs, the star must have an 
eccentricity close to unity, and it has been actively dis¬ 
cussed that the Kozai-Lidov (KL) mechanism could play 
important roles In the KL mechanism, the inner 

eccentricity of a hierarchical triple system increases due 
to the exchange of inner and outer angular momenta Q . 

In this paper, concentrating on the MBH inspirals, we 
theoretically examine the time evolution of the hierar¬ 
chical triple systems that specifically satisfy the follow¬ 
ing two conditions; (I) the KL mechanism is initially 
suppressed by the inner apsidal precessions (due to non- 
Keplerian potentials (§-[lOl|). but (2) it gradually becomes 
effective, along with the slow (to be more precise, adia¬ 
batic) contraction of the outer orbits (see also (Til. H^). 

Our primary aim is to develop a simple geometrical 
approach for the adiabatic evolution of the KL mecha¬ 
nism (see also (l3j | for similar studies mean-motion reso¬ 
nances). More specifically, we analyze the phase-space 
structure of the corresponding time-dependent Hamil¬ 
tonian, and concurrently apply the arguments of the 
adiabatic invariance. This approach would be comple¬ 
mentary to numerical studies, for understanding the KL 
mechanism of the triple systems including MBH binaries. 
The KL mechanism has been also applied to various as¬ 
trophysical phenomena [sl-ITol Hdl - I^ , and our approach 
could be useful for studying them. 

This paper is organized as follows. In section H, we an¬ 
alyze the hierarchical systems composed purely by three 
point-masses, nevertheless including relativistic effects. 


These simple systems are suitable not only for explain¬ 
ing our geometrical approach but also for numerically 
examining its validity. Then, in section HI, we apply our 
approach to individual stars in a stellar cluster associ¬ 
ated with slowly shrinking MBH binaries. We find that 
the eccentricities of the stars show intriguing evolutions. 
These results could be well explained by a peculiar bifur¬ 
cation pattern caused by the Newtonian potential of the 
stellar cluster. We also discuss the effects of relaxation 
processes for our geometrical approach. 

II. KL MECHANISM FOR PN SYSTEMS 
A. Numerical Experiments 

We first examine the evolution of relativistic hierarchi¬ 
cal triple systems, only including gravitational interac¬ 
tion with the Post-Newtonian (PN) approximation. This 
is for demonstrating our geometrical approach and also 
providing an example in distinction from the Newtonian 
effects of a stellar cluster (discussed in section HI). 

Our triple systems are made from two spinless MBHs 
(masses: toq, m 2 ) and a stellar mass object (mi). The 
inner binary is composed by mo and mi, while the re¬ 
maining MBH m 2 is the outer tertiary. Below, we use 
the geometrical unit with c = G = rriQ mi -I- m 2 = I 
and the labels j G {1,2} for the inner and outer orbital 
elements. We denote the semimajor axes by Uj, the ec¬ 
centricities by Cj, the arguments of pericenters by gj and 
the inclination between the inner and outer orbits by I. 
For the inner angular momentum, we also define 

Gi = (1) 

and 

Ji = Gi cos I (2) 

(the component normal to the outer orbital plane, Ji < 
Gi < 1). To simplify our arguments, we fix mo = 0.3, 
mi = 0.3 X 10“® and m 2 = 0.7 

To begin with, we make numerical simulations for the 
relativistic triple systems. For the direct three-body cal¬ 
culations, we use the equations of motion derived from 
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the 2.5PN ADM Hamiltonian [2l|, formally expressed as 
HadM = + HipT<j + i?2PN + Y ^^2.5PN (3) 


Here iJN,^^iPN and i? 2 PN are the Newtonian, IPN, and 
2PN terms. The last one i? 2 . 5 PN is the leading dissipative 
term and Y is an artificial parameter to speed up our 
numerical calculations. 

Using a fourth-order Runge-Kutta scheme with Y = 
104 0 , we numerically solve the evolution of the follow¬ 
ing three systems, s-I, s-II and s-IH. Their initial con¬ 
ditions are ( 614 , 02 , 1 , 62 , 1 ) = (0.001,6300,0.0001) for s- 
I, (0.001,6300,0.3) for s-fi and (0.3,5700,0.0001) for s-IH. 
Other orbital elements are commonly set at ai,i = 400, 
5 i,i = 7’'/2 and 32,1 = 0. We also choose the initial inclina¬ 
tion li to realize Ji,i = — e\^ cos I\ = 0.150 {I\ ^ 80° 

for our systems). Our main targets are s-I and s-M with 
different initial inner eccentricities 61 , 1 . The system s-ll is 
for briefly studying effects of the outer eccentricities 62 , 1 . 
Note also that we set the above orbital parameters for 
demonstration of our geometrical approach. More realis¬ 
tic systems (but more complicated for direct calculations) 
would be analyzed in the next section. 

We terminate our calculations, when the instantaneous 
separation between rriQ-mi becomes smaller than IOtoq. 
All the three systems are dynamically stable until this 
termination condition (satisfying the stability criterion 
of HI). 

Here, we summarize the characteristic timescales. For 
the outer orbits of s-I and s-IH, we have the merger time 
'72,g\^ 4.82 X I 04 O(a 2 / 6000 ) 4 ( 104 /r) due to GW emis¬ 
sion [23j. In the following, we use the decreasing dis¬ 
tance 02 as an effective time variable. Meanwhile, in the 
range 02 > 4000 relevant for the analysis below, we can 
neglect the radiation reaction for the inner orbits with 
'Ui,GW > 6072 ,GW- The timescales of the KL mechanism 
and the PN apsidal precession are also important for the 
inner binary. Their typical values are 


Tkl = 
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(4) 
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FIG. 1. Evolution of the inner orbital elements ei and g\, 
as a function of decaying 02 . Panel lA: results from the di¬ 
rect three-body calculations for s-I (red points) and s-H (green 
points). Panel IB: results from the direct three-body calcu¬ 
lation for s-IH (red points) and the corresponding analytical 
predictions (blue points) for the minimum and maximum of 
ei at each a 2 . 


evolution of Ci and gi as a function of 02 . In the case of 
s-I, 61 suddenly starts to increase at 02 ^ 5200 (showing 
the characteristic feature of a pitchfork bifurcation [^ 1 . 
and, concurrently, gi turns from circulation to libration. 

For s-IH, 61 initially oscillates, but it momently van¬ 
ishes at 02 ^ 4800. Subsequently, ei starts to increase, 
and <71 settles into libration. In Fig. lA, the evolution of 
the system s-II is similar to s-I. Therefore, the basic as¬ 
pects of the transition could be understood by analyzing 
the simpler system s-I with 62 = 0 (as we limit below). 
In common with the three systems, the mean values of 
61 evolve slowly with time. 

It has been known that, for 62 0, the inner eccentric¬ 

ity 61 can experience a resonant-like excitation around 
the epoch with Tkl Pipn SIS- Recently, Liu, Lai 
and Yuan reported that this excitation is caused by 
an apsidal precession resonance wi — W 2 — 0. Here, 
ii 7 i(= rii + gi) is the longitude of the pericenter and 
is the longitude of the ascending node. 


where Pi and P 2 are the inner and outer orbital periods 
Q. As we have P 2 ,gw Pkl,Pipn for our three sys¬ 
tems, the outer orbits shrink adiabatically. Also, the KL 
mechanism is initially suppressed by the inner relativistic 
precession with Tkl > Pipn- 

Next, we describe our numerical results for s-I and 
s-in. As expected, the shrinkages of oi are small with 
|Jai/ai,i| < 0.04 down to the epoch of 02 = 4000 (fluctu¬ 
ations of Ji: < 0.1). Figure. 1 represents the 


^ We obtained almost the same results with Y = 500, by dropping 
the time-consuming 2PN term. 


Our primary target in this paper is the inner orbital 
evolution for 62 = 0. But, here, we briefly comment on 
this resonance. For the run s-H mentioned above, we 
have dji ^ 1772 due to the hierarchies of the masses and 
orbits. Therefore, the resonant state tui — 1372 ~ 0 was 
not realized in the run. 

In the next section, we study a more realistic system 
including a stellar cluster potential with 62 = 0. It would 
be Interesting to examine the effects of the apsidal preces¬ 
sion resonance for 62 0. But, such an extension is far 

beyond the scope of this paper, since we need additional 
inputs, e.g. a detailed modeling of the outer precession 
rate fi72. We remain this issue as a future work. 
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B. Averaged Hamiltonian 


C. Geometrical approach 


To geometrically explain the interesting behaviors in 
Fig. 1, we briefly describe the employed Hamilto¬ 
nian. For studying the long-term evolution of the triples, 
we take the double averages of the original three-body 
Hamiltonian with the inner and outer mean anomalies 
After approp riate scaling, we obtain the relevant 
Hamiltonian |lfll. 


'Ht, 1PN = Hqp -b "Hipn (6) 

with 

Hqp = -3G? - 15^ - 15(1 - Gl) (^1 - cos2gi, 

(7) 

7^ipn = -^- (8) 

Here, Hqp is the Newtonian quadrupole term due to the 
outer tertiary, and is obtained by perturbatively evalu¬ 
ating the tidal field of the tertiary m 2 with the expan¬ 
sion parameter ( 01 / 02 ) <C 1. The subsequent-order (oc- 
tupole) term vanishes for 62 = 0. The IPN term "Hipn 
leads to an apsidal precession of the inner orbit [T^ . 

The dynamical degree of freedom of this Hamiltonian 
is only the canonical variables {gi,Gi). Below, as sup¬ 
ported by the numerical experiments, we put oi = const 
and Ji = const (ignoring the inner radiation reaction) 0 . 
The decreasing parameter 7 is defined as 


16(mo -I- mi )^02 ^ / Tkl \ 

(m 2 o|) of Tipn/ ’ 


(9) 


and contains the information of the outer orbits in 
"Ht.ipn- The IPN effect gradually becomes weaker, rel¬ 
ative to the tidal field of the infalling tertiary. 

In the (gi, Gl) coordinate, the angle variable gi is sin¬ 
gular at Gl = Ji and 1, since gi loses its geometrical 
meanings there (either I = 0 or ei =0). These coordi¬ 
nate singularities should be regarded as two degenerate 
points, and can be appropriately handled with the canon¬ 
ical transformations Q; 

{x',y') = V 2 (Gi - Ji)(cos 5 i,singi) ( 10 ) 
{x,y) = v^2(l - Gi)(cosgi,singi) (11) 

As expected from geometric symmetry, these two points 
are fixed points in the phase space, irrespective of the 
parameter 7 . 


^ Also, our Hamiltonian 'Kt.ipn is independent of their conjugate 
variables. 


For our geometrical approach, we need to understand 
how the phase-space structure changes along with the 
slowly decreasing parameter 7 . From general arguments 
on Hamiltonian systems, a fixed point must be either a 
center (stable fixed point) or saddle (unstable fixed point) 

H- 

In Fig. 2, we show the phase-space structure of the 
Hamiltonian 'Ht,ipn( 3 i, Gi; 7 ) for 7 = 14.89,10.03 and 
7.37. First, we focus on the solid black curves (con¬ 
tour samples of "Ht.ipn), the black points (unstable fixed 
points), and the red broken curves (separatrixes). These 
elements evolve in the following order. 

(PI) For 7 > 12 — 20J/ = 7 c (Fig. 2A), we only have 
circulating trajectories. Accordingly, the fixed points at 
Gl = 1 and Ji are both stable. The latter remains stable 
during whole. 

(P2) At the incident 7 = 7 c (between Figs. 2A and 
2 B), the fixed point Gi = 1 transits from a stable point 
to an unstable point (pitchfork bifurcation M)- For 
Ji = 0.150, this bifurcation is at 7 c = 11.55 (nicely cor¬ 
responding to 02 — 5238 in Fig. lA). 

(P3) For 7 < 7 c (Figs. 2B and 2C), we simultaneously 
have circulating and librating trajectories, segregated by 
the separatrixes. As 7 decreases, the stable points at <71 = 
± 7 r /2 take smaller Gi. The area inside the separatrixes 
increases monotonically. 

Then, how can we trace the evolution of individual 
systems in the sequence of Figs.2A-2C? Here, we explain 
a powerful tool. Let us consider a periodic motion de¬ 
scribed by a Hamiltonian including a param¬ 

eter A. If the parameter A changes more slowly than the 
period of the motion, the area surrounded by the periodic 
trajectory S = <f pdq is conserved and called an adiabatic 
invariant [261127|| . 

Now, let us trace the evolution of the system s-III in the 
phase-space, based on the adiabatic invariant. As in Fig. 
2A, this system initially circulates in the ( 51 , Gi) plane, 
but it is periodic in the {x, y) plane with the calculated 
adiabatic invariant S = 0.129 = S'!!!. Here the area Sni 
corresponds to the blue region in Fig. 2A. Because of 
the area conservation, at the stage of Fig. 2B, s-M is 
expected to stay circulating on the thick blue curve. 

But, as 7 decreases, the total area inside the two sepa¬ 
ratrixes increases, and it eventually becomes identical to 
S'ni between Fig. 2B and 2C (around 7 ^ 8.89). At the 
separatrix crossing, ei momently vanishes, and the sys¬ 
tem s-in starts to librate either around gi = 7r/2 or — 7 r /2 
(equal probability from the symmetry), due to the conti¬ 
nuity of the trajectory. Associated with this disjunction, 
the adiabatic invariant S falls to one-half iS'n/2. But this 
gap is not contracting with the adiabatic invariance, since 
the period becomes infinite for the trajectory just on the 
separatrix. After crossing the separatrix, we can again 
follow the system, using the area conservation around 
one of the two stable points. For example, the averaged 
value of ei increases with time. In this manner, we can 
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FIG. 2. PN evolution of the phase-space of the inner binary with Ji = 0.150 (plotted in the limited range 0.75 < Gi < 1). 
Figs. 2A, 2B and 2C correspond to 7=14.88, 10.03 and 7.37 (decreased by the tertiary). The following two types of information 
are presented simultaneously in each panel. [1. the overall phase-space structure of the Hamiltonian Ht,ipn] solid black curves: 
contour samples of Ht.ipn, black points: stable fixed points, red dashed curves: separatrixes, and arrows: directions of flows. 
[2. evolution of the specific system s-III] The blue regions have equal area ^ni = 0.129 (corresponding to the adiabatic invariant). 
The triangles represent the results from the direct three-body calculation for s-IH. 


geometrically understand the evolution of s-DI. 

Until now, we have discussed the phase-space mapping 
of the specific system s-III. Generally, in the mapping 
from Fig.2A to 2C, a system initially closer to Gi = 1 is 
transported to a librating trajectory closer to one of the 
two stable fixed points (crossing the separatrix earlier). 
As the present mapping is thus simple and symmetric, 
we can justify the area conservation also by applying Li- 
ouville’s theorem to certain group of systems {e.g. these 
composing the blue region in Fig. 2A). Indeed, in Fig. 2, 
the total area of the blue region is the same, before and 
after the separatrix crossing. 

Next, we examine the validity of our approach, using 
the direct three-body calculations. In Fig. 2, the trian¬ 
gles represent the numerical results for s-DI. The associ¬ 
ated area is conserved at high precision. Assuming the 
area conservation, we can inversely predict the evolution 
of ei easily. In Fig. IB, with the blue points, we repre¬ 
sent the maximum and minimum values of ei at each 02 . 
They well reproduce the direct calculations 0. 

In the case of s-I (with much smaller area Si <C -S'n), 
the trajectory turns into libration soon after (P2) and 
it stays close to the stable fixed point throughout Fig. 
2B and 2C. This is why the distinctive features of the 
pitchfork bifurcation (originally associated with the fixed 
points) appear in Fig. lA. 


with an outer infalling tertiary MBH m 2 . We deal with 
component stars for which (i) the dominant force is the 
Newtonian attraction by the central MBH mo, but (ii) 
its IPN correction would be smaller than the Newtonian 
potential of the cluster itself. Here, we ignore the IPN 
effect for simplicity, but it could play important roles 
for stars with small pericenter distances (see e.g. i). 
In this section, we consider outer MBH that have larger 
separation (still on circular orbit) than the previous PN 
systems (see §HI.B for parameters in the physical units). 
We do not specify the physical processes for the outer 
orbital decay, but the infall rate is assumed to be suf¬ 
ficiently small {e.g. after the stall of dynamical friction 
i). We also disregard dissipative effects for the inner 
orbit {i.e. setting oi = const and Ji = const). The re- 
laxation processes would be discussed in §111.B (see e.g. 

smiig). 

A. Geometrical approach 

We first present the averaged Hamiltonian. The cluster 
is initially assumed to be spherical with the local density 
profile p{r) = pi{r/ai)~^ (/3 < 3) around the distance 
r = ai in interest (pi: the mean density there). Ignoring 
the back-reaction (time variation) to the stellar potential, 
the scaled Hamiltonian is written as 


III. NEWTONIAN POTENTIAL OF STELLAR 
CLUSTERS 

We next apply our approach to the evolution of indi¬ 
vidual stars in a stellar cluster around a central MBH mo 


'Ht.sp ='Hqp-b "Hsp- (12) 


Here, 'Hsp is the term due to the stellar potential 


"Hsp = 


I + /3(-l-f/3) 



(13) 


® At the very end of the run, there exist mismatches caused by 
breaks-down of the averaging method for ei 1 


which is obtained by perturbatively expanding a hyperge¬ 
ometric function with ei = i/l — Gf [ 3 . In our geomet¬ 
ric approach, we neglect the small higher-order correction 
o(ef) that vanishes identically at /3 = 0 and 1 (though 

















5 


exactly handled later with /3 = 3/2). Namely, we put 
"Hsp = ^(1 ~ Gi) in terms of the canonical variable Gi. 
Then, all the information of the cluster is included in the 
single decreasing parameter 


V = 


167rpia2 

m2 


(14) 


with 0,2 = 9(02 < 0 . 

The precession timescale due to the stellar potential 
TsP i is 


Tsp 


(3 - 9)mo 
iwpiaf ^ 


(15) 


Then, we can represent rj by using Tkl and Tgp as 


127r Tkl 

3 - /3 Tsp ■ 


(16) 


Next, we explain the evolution of the phase-space 
structure, following Fig. 3 (see also i)- The most re¬ 
markable point below is the bifurcation (P4’) at r]c 2 = 12 
(independent of Ji!). 

(PI’) For rj > —18 -I- 30/J/ = rjci (Fig- 3A), we only 
have circulating trajectories in the retrograde direction. 

(P2’) At T] = rjci (between Figs. 3A and 3B), we have 
a pitchfork bifurcation similar to (P2). The fixed point 
Gi = Ji turns from a stable point to an unstable point ( 
in contrast to Gi = 1 for (P2), see also Appendix A). The 
new stable point (black point in Fig. 3B) is generated at 
9i = 7r/2. 

(P3’) For 77 c 2 = 12 < ?7 < ryci (Fig. 3B), we have 
librating trajectories inside the separatrix (red dashed 
curve). The inner area increases with time. 

(P4’) At T] = r]c 2 (Fig. 3C), due to a cancellation 
between Hqp and Hsp, the Hamiltonian "Ht.sp becomes 
independent of Gi on the two vertical lines gi = 0 and tt. 
This is a quite peculiar state in which fixed points spread 
fully on the two lines. At this moment, the stability of the 
two points Gi = Ji and 1 interchange, and only librating 
trajectories exist. 

(P5’) For 7 ?e 3 = -18 + 30J2 < g < g ,2 (Fig. 3D), 
we again have both librating and circulating trajectories. 
The latter now have prograde motion. The area inside 
the separatri x dec reases. 

For Ji > y^3/5, we have r?c 3 > 0. Thus, if g decreases 
down to ryc 3 , the stable point at gi = 7r/2 disappears, 
and the fixed point Gi = 1 turns to be stable (inverse of 
the bifurcation at (P2)). We then have only circulating 
trajectories. 

Now, using the area conservation, we follow the evo¬ 
lution of two representative stars s-IV and s-V in the 
time sequence of Fig. 3. Their initial conditions ( 51 , 
Gi) at 77 = 200 are (0,0.95) for s-IV, and (0,0.6) for s- 
V both with Ji = 0.5. In Fig. 3A, we can effectively 
regard the areas of the colored regions (blue: s-IV and 
green: s-V) as conserved quantities. Conserving these 
areas, the two systems would evolve as shown with the 
thick blue/green curves in Figs. 3B-3D. Note that, at 


the first separatrix crossing during (P3’), the maximum 
eccentricity (i.e. Gi = Ji) is surely realized (in contrast 
to Gi = 1 for the IPN case). 

Interestingly, the magnitude of their eccentricities is 
inverted between Fig. 3A and 3D. This is because of the 
passage of the peculiar state (P4’). Consequently, nearly 
circular orbits (Gi ^ 1 ) would be quickly converted to 
eccentric ones (Gi ^ Ji) around 77 = 12. This sharp tran¬ 
sition is quite different from the more gradual evolution 
in the PN case, and would be less likely to be prevented 
by other processes. 

If the cluster initially has an isotropic velocity distribu¬ 
tion, the stars accordingly have a homogeneous density 
distribution in the ( 51 , Gi) plane [3l[. Then, because of 
Liouville’s theorem, the homogeneity is unchanged in the 
subsequent stages in Fig. 3, and our assumption about 
the stationarity of the cluster potential would be a good 
approximation. Meanwhile, for an initial distribution bi¬ 
ased to small eccentricities, the secure realization of the 
maximum eccentricity during (P3’) could help to refill the 
loss cone, and might be very important astrophysically. 

In Fig. 3, we additionally present the results obtained 
by numerically solving the Hamilton equations; 

^ 97It,sp . 9'Ht,sp (1-7^ 

with (3 = 3/2 and g = —10“^ (including o(ef)-term for 
"Hsp)- Except for some fluctuations around ry = 12, our 
simple geometrical predictions fit the numerical results 
well. 


B. Effects of Relaxation 

Until now, we have neglected the diffusion of the energy 
and angular momentum associated with the discreteness 
of the stellar distribution. In this subsection, we estimate 
the characteristic timescales of the relevant relaxation 
processes. 


1. Model parameters 

To begin with, we set the fiducial model for the overall 
cluster density profile as 


p(r) 


3too / 
dTrrj^ V^h/ 


(18) 


with the slope /3 = 3/2 [32j. Here is the influenced 
radius = Gm^la'^ with the one-dimensional velocity 
disp ersion a that is assumed to satisfy the Af, - a relation 


M.(= mo) = 2.09 X 10® 


200 km/s 


5.64 

Mq. 


(19) 
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FIG. 3. Evolution of the phase-space for stars in a stellar cluster around a MBH (hxed at Ji = 0.50). Figures 3A, 3B, 3C 
and 3D are for 7^=200, 20, 12 and 5 decreased by an infalling tertiary MBH. Because of the symmetry, we only plot the range 
0 < (/I < TT. As in Fig. 2, the following two types of information are simultaneously plotted. [1. the phase-space structure of 
'Ht.sp] The symbols (solid black curves, black dots, red dashed curves and arrows) are defined in the same manner as Fig. 2. 
[2. evolution of the two representative stars s-IV (blue) and s-V (green)]. The areas of the colored regions are respectively the 
same. The triangles and squares are obtained by numerically integrating the Hamilton equations with /3 = 3/2 and »)=—10~®. 


From Eqs. m and (HU), the timescale Tsp of the mass 
precession is given by 


without depending on oi (for /3 = 3/2). 

In the standard picture of MBH binary evolution Q, 
the orbital decay timescale Tuec = 02 / 0-2 is considered to 
increase sharply around the hard-binary separation 


Oh = 


m2 rh 


mo - 
= 3.50 


m 2 4 

’ (- 


niQ 


g -1- 1 V IO8M0 


0.645 


pc 


( 21 ) 


with q = m 2 /mo < 1. The subsequent orbital decay 
would strongly depend on the exterior of the MBH binary 
{e.g. triaxial potential [13,H^), and the actual timescale 
Tdec is still under active discussion. We are specifically 
interested in the orbits of stars around the peculiar bifur¬ 
cation at ?7 = 12. For applying the adiabatic invariance, 
the decay timescale T^ec must be smaller than the relax¬ 
ation timescales, but larger than Tkl that is explicitly 
given as 


Tkl 


5.45 X 10® / mo \ 
~q VlO^Mej 


X 




( 22 ) 


2. Resonant Relaxations 

Now, we discuss the timescales for the relaxation pro¬ 
cesses. As is well-known for nuclear star clusters, the 
timescale for the two-body relaxation is generally much 
larger than that for the resonant relaxation. Therefore, 


in the following, we concentrate on the resonant relax¬ 
ation. 

We firstly summarize the basic aspects of the resonant 
relaxation for a conventional system composed only by a 
central MBH and the associated stellar clusters Q . Then, 
we discuss how the resonant relaxation would be modified 
by the outer MBH m 2 . 

In a spherical stellar cluster around a central MBH 
mo, the individual stars maintain their Keplerian orbits 
for the precession timescale Tsp. The scalar resonant 
relaxation is the long-term accumulation of the gravita¬ 
tional interaction between these Keplerian orbits of the 
short-coherence time Tsp. Both the magnitude and ori¬ 
entation of the angular momenta of the stars change with 
the characteristic timescale Q 




mi 


= 1.6 X 10® 



mp 

IO^Mq 


1/2 


Ql \ 

0.3pc J 


3/2 


yr- 


(23) 


Meanwhile, the vector resonant relaxation should be 
regarded as interaction between the mass annulus ob¬ 
tained after the orbital averaging over the precession 
timescale Tsp. It changes only the direction of the angu¬ 
lar momentum (due to the symmetry) with the charac¬ 
teristic time 0 


Trr.s 

y/]v’ 


(24) 


where N is the number of stars around the semimajor 
axis in interest. 

Next, we discuss the impacts of the additional outer 
MBH m 2 on the two-types of resonant relaxations. Now, 
the individual Keplerian orbits maintain their shapes for 
min[TKL,T sp]. Therefore, for Tkl ^ Tsp, the timescale 
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of the scalar resonant relaxation would be a factor of 
Tsp/^kl = (’7/12)”^ tt/ (3 — /3) times larger than in 
Ea. (l23l) . This factor is of order unity for 77 ^ 12, 

Here, it is important to note that the outer MBH also 
generates precession of the longitude of the ascending 
node Hi of the stars with the characteristic timescale Tkl 
(because |Hi| = |9 Ht,sp/9Ji| -- |97^t,sp/9Gi| = |gi|). 
Therefore, around 77 ~ 12, the vector resonant relax¬ 
ation should be regarded as gravitational interaction be¬ 
tween axisymmetric (torus-like) mass distributions with 
the common symmetric axis (rather than the inclined an¬ 
nulus for the conventional case). But there is no trans¬ 
fer of torque between them (due to the symmetry), and, 
consequently, the vector resonant relaxation becomes in¬ 
effective around the regime 77 ^ 12 . 

We conclude that, even though the decay timescale 
Tdec = 02/02 is uncertain at 02 < Oh, it should be smaller 
than Trr s in Eg. (1231) for applying the adiabatic invariance 
to the bifurcation around 77 = 12 . 

IV. DISCUSSIONS 

We have developed a geometrical approach for the adi¬ 
abatic evolution of hierarchical triple systems including 
MBH binaries on circular orbits. We analyzed the basic 
profile {e.g. fixed points and separatrixes ) of the phase- 


space structure of the corresponding Hamiltonian. In 
addition, we apply the arguments of the adiabatic invari¬ 
ance to follow the time evolution of individual systems. 

At first, we demonstrated the validity of our approach 
using a simple relativistic systems that can be also exam¬ 
ined by direct numerical calculations. Then, we analyzed 
more realistic systems, star clusters associated with inspi- 
raling MBH binaries. We found that the eccentricities of 
the stars could evolve quite interestingly. This is caused 
by the peculiar bifurcation patterns induced by the New¬ 
tonian potential of stellar clusters, in sharp contrast to 
the precedent relativistic systems. Our simple geomet¬ 
rical approach enables us to clearly understand how the 
inner orbits are deformed, in response to the time vari¬ 
ation of the related phase-space structure. It would be 
applicable to the KL mechanism in various astrophysical 
contexts, possibly including viscous effects to the inner 
orbit at least for a short time. 
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equivalent to 51 = 0 for the apsidal precession rate, rep¬ 
resenting the balance between the quadrupole effect and 
the competing IPN or SP effect. In concrete terms, we 
collectively have 

Q(Gi) = -36Gi + 60^ = -fcGr“. (A2) 


Appendix A: Generation of a New Fixed Point 


As shown in Figs. 2 and 3, the phase-space evolution 
is largely different between the two Hamiltonians Ht.ipn 
and 'Ht,sp given in Eqs. ( 6 ) and ( 12 ). When relatively 
increasing the quadrupole term (ie. ^,r] = -|-oo —>■ 0 ), 
one of the remarkable contrasts is the emergence of the 
new fixed point either at Gi = 1 (the epoch (P2) for 
IPN) or at Gi = Ji (the epoch (P2’) for SP) both with 
gi = ±.7:12. In this appendix, we discuss how this con¬ 
trast is related to the sign of the apsidal precession (IPN: 
prograde and SP: retrograde). 

In the present case, the condition for the fixed point is 
formally expressed as 




dGi 


gi=Trj2 


= 0 


(Al) 


for 'Ht' = 'Ht.ipn and 'Ht,sp- This equation is physically 


Here, the function Q{Gi) originates from the quadrupole 
term and monotonically decreases in the relevant range 
Ji < Gi < 1. The right hand side shows the apsidal 
precession rates induced by the IPN (a = 2 and k{= 
37 ) > 0) and the SP (a = — 1 and k{= —2rf) < 0). 

Now, for simplicity, we assume Ji < y^3/5 (identical to 
the condition for the standard KL-libration). Then, we 
have Q{Gi = Ji) > 0 and Q{Gi = 1) < 0. Furthermore, 
with respect to the two types of variations k = -l-oo —>■ 0 
or fc = —00 —>■ 0 , either of two solutions; Gi = Ji or 
Gi = 1 initially appears in Ea. (IA2l) . irrespective of the 
index a G [—1,2]. Indeed, depending only on the sign of 
fc, we have Gi = 1 (for fc = +00 —>■ 0 ) and Gi = Ji (for 
fc = —00 —>■ 0). 

This result explains the sharp contrast between (P2) 
and (P2’). Here, the sign of the apsidal precession is 
crucially important. Note also that, even if a small cor¬ 
rection term is added to the right hand side of Eq. (IA2I) , 
this robust result is unchanged. 




